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SUMMARY 

Numerically solving the incompressible Navier-Stokes equations is known to be time- 
consuming and expensive. Testing of the INS3D computer code, which solves these equations 
with the use of the pseudocompressibility method, shows this method to be an efficient way to 
obtain the steady-state solution. The effects of the waves introduced by the pseudocompressibility 
method are analyzed and criteria are set and tested for the choice of the pseudocompressibility 
parameter which governs the artificial sound speed. The code is tested using laminar flow over a 
two-dimensional backward-facing step, and laminar flow over a two-dimensional circular cylinder. 
The results of the computations over the backward-facing step are in excellent agreement with 
experimental results. The transient solution of the flow over the cylinder impulsively started 
from rest is in good agreement with experimental results. However, the computed frequency of 
periodic shedding of vortices behind the cylinder is not in agreement with the experimental value. 
For a three-dimensional test case, computations were conducted for a cylinder-end wall junction. 
The saddle point separation and horseshoe vortex system appear in the computed flow field. The 
solution also shows secondary vortex filaments which wrap around the cylinder and spiral up in 
the wake. 


INTRODUCTION 

Computational fluid dynamics (CFD) is rapidly advancing both as a tool for theoretical 
study and as a tool for engineering design. This is partly due to the increasing power and 
efficiency of the high-speed supercomputer. It has, therefore, become feasible to numerically 
simulate complicated fluid-flow phenomena associated with complex geometries. Thus, one of the 
primary pacing items in computational fluid dynamics research is the development of new solution 
methodologies (ref. 1). This paper is the result of research at NASA Ames Research Center in 
coordination with the development of a three-dimensional, incompressible, Navier-Stokes solver 
utilizing primitive variables and generalized curvilinear coordinates. 

PURPOSE OF THE STUDY 

There are many important incompressible flow problems in engineering applications, 
some of which include hydrodynamics, liquid fuel flows in propulsion, and low-speed air flow. 
When the flow is two-dimensional in nature, there are several efficient methods for numerically 
solving the incompressible Navier-Stokes equations. These flow solvers do not signifi r antly impact 
a supercomputer’s memory and time constraints because two-dimensional flow fields have small 
data base requirements, e.g., the stream function-vorticity formulation is one popular method 
(refs. 2-5). However, many incompressible viscous flows are highly three-dimensional in nature 
and there is no staightforward three-dimensional extension to this method; therefore, the use 
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of primitive variables (velocities and pressure) is preferred. There are several three-dimensional 
codes which solve the compressible Navier-Stokes equations, examples of which are given by 
references 6-9. Since the speed of sound approaches infinity at the incompressible limit, it is 
not computationally efficient to use these codes for computing incompressible flows. Therefore, 
there is a need for an efficient three-dimensional incompressible Navier-Stokes solver which can 
be used for a realistic three-dimensional geometry. 


SOLVING THE INCOMPRESSIBLE NAVIER-STOKES EQUATIONS 

One of the major concerns in numerically solving the incompressible Navier-Stokes equa- 
tions in primitive variable form is the method used for evaluating the pressure. The pressure does 
not explicitly appear in the continuity equation yet it acts as a dynamic parameter in ensur- 
ing that continuity is satisfied, and that the velocity field is divergence-free. One method that 
was developed by Harlow and Welch (ref. 10), for use with their MAC (marker and cell) code, 
involves the solution of a Poisson equation for pressure. By taking the divergence of the mo- 
mentum equations, a Poisson equation in pressure is derived. In using this method, the normal 
time-advancement would be as follows: the momentum equations are solved for the velocities, 
then the Poisson equation is solved for the pressure which ensures a divergence-free flow at the 
next time-step. The solution of the Poisson equation usually requires a relaxation scheme to 
iterate on pressure until the divergence-free condition is reasonably satisfied. This method has 
been widely used, particularly for solving two-dimensional problems, (refs. 11 - 17). However, 
the solution of the Poisson equation can be prohibitively time-consuming for three-dimensional 
problems. 

In his 1968 paper, Chorin (ref. 18) presents the basis for applying the fractional step 
method in solving the incompressible Navier-Stokes equations. The theoretical basis can be found 
in another paper by Chorin (ref. 19) and in the book by Temam (ref. 20). The essential mechanics 
are as follows: in the computation the time-step is split into two levels. In the first level, a modified 
set of momentum equations are solved for an auxiliary velocity field. The modification to the 
momentum equations is made by removing the pressure gradient. All that remains is to solve for 
the pressure which will map the auxiliary velocity into a divergence-free field. For this Chorin 
produces two equations, one in which the new velocity is calculated from the auxiliary velocity 
minus the pressure gradient, and in the second equation the change of pressure is calculated from 
the divergence of the new velocity. These have the dimensionless form 


= A<^p" +1 

(1.1) 

= p n — ADiv(u n+1 ) 

(1.2) 


where A is a constant chosen to facilitate convergence. These can be solved by iteration and result 
in a divergence-free flow field. By substituting equation (1.1) into equation (1.2) it is seen that 
this is equivalent to solving the Poisson equation in pressure 


V*p=^Div(„““') (1.3) 

If the convective and viscous terms from which the auxiliary field was computed were substituted 
into equation (1.3) then the resulting Poisson equation would be that of the previous method. So it 
can be seen that by performing the decomposition of the momentum equations and by computing 
the auxiliary field first results in a savings in the required number of floating point operations in 
each iteration. This method has been used by Kim and Moin (ref. 21) with some success to write 
a three-dimensional code in Cartesian coordinates which solves the Poisson equation in equation 
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(1.3) with a direct method. However, it is still preferable to avoid solving a Poisson equation in 
three-dimensional curvilinear coordinates as this would take more computing time than is desired. 

To avoid the solution of a Poisson equation in pressure, there is the artificial compress- 
ibility method as first developed by Chorin (ref. 22). In this method, a time derivative of pressure 
is added to the continuity equation 


dp 

dt 



= 0 


(1.4) 


Together with the momentum equations this forms a hyperbolic-type system of equations. The 
equations can then be simultaneously advanced in time using a suitable algorithm. As the so- 
lution converges, the time derivative of pressure approaches zero and continuity is satisfied in 
the steady state. The artificial compressibility introduces pressure waves of finite speed into the 
incompressible flow, whose speed of sound would otherwise be infinite. These waves die out as 
the steady state is reached. This method was combined by Steger and Kulter (ref. 23) with an 
implicit approximate factorization scheme of Beam and Warming (ref. 24) for the computation 
of vortex wakes. This combination is known as the method of pseudocompressibility and is found 
to be computationally efficient. Because the continuity equation is not satisfied until the solution 
is converged to a steady state, this method is not time-accurate. But by temporarily relaxing 
the requirement of continuity, the steady state solution can be obtained at a significant sa\ ings 
in computing time. The use of this method currently appears to be the most viable approach to 
obtain incompressible steady state solutions for a complex three-dimensional geometry. 

Pseudocompressibility was chosen by Kwak et al (refs. 25-26) in developing an efficient 
and robust code called INS3D, a three-dimensional, incompressible viscous flow solver in gener- 
alized curvilinear coordinates. This code was the main tool used in this research in an effort to 
study the best method for solving the incompressible Navier-Stokes equations. 

The outline of this papers is given below. In the second section, the governing equations 
will be presented. The finite-difference algorithm will be described in the third section and the 
method of pseudocompressibility will be discussed further in the fourth section. Computed results 
will be presented in the fifth section for several test problems. These include a two-dimensional 
backward-facing step, flow behind a two-dimensional circular cylinder, and the flow around a 
cylinder-plate junction. 


GOVERNING EQUATIONS 


Incompressible Navier-Stokes Equations 

The equations governing unsteady, three-dimensional, incompressible flow with constant 
density are presented. In these equations t is time; x, y, and z are the Cartesian coordinates, u, v, 
and w are the corresponding velocity components; p is the pressure; and r t , is the viscous s ress 
tensor. The equations have been written in dimensionless form using the following dimensionless 
variables 


u = 
x = 
t = 


u 

V XV 
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The subscript ref denotes reference quantities and in the following equations the tildes (') have 
been dropped for convenience. 

The momentum equations are given by 


dt " dx 

and the continuity equation is 
where 


d A d . , d . . d 

q + ~ €v ) + dy( f ~ + dz^~ &) = 0 


du dv dw 
dx + dy + dz ^ 


(2.2) 


(2.3) 


9 = 


u 

V 

w 



u 2 + p 


• 

vu 


wu 

e ~ 

uv 

f = 

V 2 + p 

9 = 

wv 


UU) 


vw 


w 2 + p 


(2.4) 


- 


- - 


P - 

T XX 


T yx 


^ ZT 

T xy 

fv = 

T yy 

II 

T *y 

m ^ xz _ 


I 

N 


- TzZ _ 


To implement the pseudocompressibility scheme the continuity equation is replaced with 


dp du dv dw. 

a7 + /?( ai + ^ + aI> 


0 


(2.5) 


Here 0 is the pseudo compressibility parameter which is chosen to facilitate fast convergence and 
will be discussed in more detail later. 

Introducing Cartesian tensor notation, the viscous stress tensor can be written as 

(2.6) 


2 uS 


v 


R 


tJ . 


Here, u is the kinematic viscosity, S',, is the strain rate tensor, and i? t , is the Reynolds stress. 
These are given by 

_ 1 , du t du , , 

? -> = s(3I 1 + -3^) (2.7) 


2 v dxj dx t 
zR-kk^ij 2u t Sij. 


R\) — „Rkk^ij (^-®) 

where it, = u, v, or w for t = 1, 2, or 3, respectively, and x, = x,y, or z for t = 1, 2, or 3, 
respectively. The normal component of the Reynolds stress is Rkk and i/ t is the turbulent eddy 
viscosity. By substituting equations (2.7) and (2.8) into (2.6) and including the normal stress, 
Rkk, in the pressure, the viscous stress tensor becomes 


, w du{ du,. 

T '> = {l/ + + aff 

In the remainder of this paper the total viscosity, (u + u t ), will be represented simply by v. 


(2.9) 
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By combining equations (2.2) and (2.5), the governing equations are written as one 
vector equation 

%-D+ £-{E - E v ) + £-{F - F v ) + |-( G - G v ) = 0 (2.10) 

at ox ay dz 


where 


D = 


P 

u 

v 

Lit; j 
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( 2 . 11 ) 
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Coordinate Transformation 


To perform calculations on three-dimensional arbitrarily shaped geometries, the follow- 
ing generalized independent variables are introduced which transform the physical coordinates 
into general curvilinear coordinates 

T — t 


t = £(*,!/, *,0 
V = »i(x,y,z,0 
C = ?(x,y,z,0 


( 2 . 12 ) 


The spatial and time derivatives become 


d _ d£ d dr) d_ df d d^ d 

dxi dx{ d£ dii dr) dxi d$ dii d£j 

d d d£j d 

dt dr dt d£j 


Here the Jacobian of the transformation is defined 

£z £y £z 

Vx V y r) z 
Ci Cy C z 


J = det 


d(£,h, C) 
d(x,y,z) 


(2.13) 


(2.14) 
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where 


* _ K _ dr, 
iz dx ,Vy ~dy' 

Before the transformation is applied, equation (2.10) is written 


± D + ±H, = 0; H t = 


E-E v 

F-F v 

G-G v 


(2.15) 


Applying the transformation to this equation yields 


d _ dF d dij d „ 

ai D+ ~MdT 1 D+ ar,dc i H ' = 0 


(2.16) 


After dividing through by the Jacobian, this can be decomposed into a conservative term and a 
remainder 

dr ' J ' d(, J at J dx, *' 


dr [ J ) + dtjJ dt 'dtjhdxJ 


(2.17) 


It can be shown that the nonconservative remainder term on the right-hand side of equation 
(2.17) is equal to zero. Thus the governing equation, in an expanded form is 


~(j ) + j\(,D + {,(£ - E.) + f,(f - F.) + UG - G.))} 

d 1 

+ q~{ jlVtD + Vx(E - E v ) + r) y (F - F v ) + rf z (G - G„)]} 
+ + $x{E - E v ) + $ y (F - F v ) + Cz(G - G„)]} 


(2.18) 


This can be rewritten in the form 


where 




D 


D 

7 

1 


E = 7 [ZtD + ( X E + i y F + £ Z G 


1 

J 


£tP + + Z y (3v + £ z /3w 

+ Cx(u 2 + p) + £ y vu + £ z wu 

£tV + i x UV + £y(t> 2 + p) + tzW 
L&u; + t X.UVJ + £ y vtv + £ z (u > 2 + p) J 


(2.19) 


I 



Now, the contravariant velocities, U, V, and W without metric normalization are defined as 

U = £t d" •+■ £yV + £zW 

V = r) t + r/ x u -I- + rj z w (2.20) 

W = + Cz«>. 


A A A 

Using these, E 1 , F, and G become 


e =i 


i 
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0 U + tt{p - 0 ) 

uU + i x p 

VU + £yP 

wU + £ z p 

0 V + r,t(p- 0 ) 

uV +T] x p 

vV + rjyp 
wV + r) z p 

pw + u(p-P) 

uW + f x p 

VW + $yP 

wW + q z p 


The viscous terms are 


E v = — d" £yE v H" 

F v — — [f7zi?u + *1yF v “H *?z^v] 

G v = — [fz-^u + SyF u + Cz^»v] 

By using the values for E v in equation (2.11), E v becomes 

0 


*« = j 


£x^XZ “H £yTyx “1“ £z^zx 
Wx^xy + r )y T yy + Vz T zy 

L f z rx2 + f y r yz + f z r zz J 


(2.21) 


( 2 . 22 ) 


(2.23) 


Recalling equation (2.9) and applying the coordinate transformation to it, the shear stress tensor 


is 


Ti, 


( d£ k duj 
"{dXidtk 


d£k duj \ 
+ dxjdik)' 


(2.24) 
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Substituting equation (2.24) into the equation for E v in (2.22) yields 


Ey — 
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du_ 

dx { dx t d£ k 
dr) dj, ± dv 
dx, Wx, d( k 

d±_dikdw 

dx t dx t d( k 
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dt 
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du dv dw v 

1 *• at at) 


d{ t 


d{ t ”dti 


0 
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where 


Im 


0 0 0 0 * 

0 10 0 

0 0 10 

.0 0 0 1 . 


If an orthogonal coordinate system is assumed, then 


■ V£y = 0; for i ^ j 


Then all the viscous terms become 
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(2.25) 


(2.26) 


For constant u (i.e., laminar flow), the contribution by the second term in parathensis 
in the above equation will approach zero as the steady state is reached and equation (2.3) is 
numerically satisfied. 
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NUMERICAL ALGORITHM 


Time- Advancing 

The numerical algorithm used to advance equation (2.19) in time is an implicit, nonit- 
erative, approximately factored, finite-difference scheme by Beam and Warming (ref. 24). The 
time-differencing used by this scheme is generally known as the trapezoidal rule and is given by 


pn+l 



+ 0(Ar 3 ) 


(3.1) 


where the superscript n refers to the n th time-step. The finite-difference form of equation (2.19) 

6 r D + S i {i!-i v ) + 6 n {P-F v ) + 6 s {G-G v )=0 (3.2) 

Q AAA 

where 6$ is the finite difference operator for for example. E, F, and G are given by equation 
(2.21) and E v , A,, and G v are given by equation (2.26). By substituting equation (3.2) into 
equation (3.1) and using D = DJ , one obtains 


D n+1 + _ E v ) n+l + 6r,(F - F v ) n+1 + 6<{G - G v ) n+1 ] 

2 

= ^J^'ME-EX+^F-Kr + i^G-G,)"} 


(3.3) 


The problem is to solve for D n+1 , and this is nonlinear in nature since E n+X = E(D n ~ fl ) is a 
nonlinear function of D n+l as are F n+1 and G n+1 . The following linearization procedure is used. 
A local Taylor expansion about u n yields 


E n+ 1 = E n + A n (D n+l - D n ) -I- 0(Ar 2 ) 
F n+l = F n + B n (D n+l - D n ) + 0(At 2 ) 
G n+l = G n + C n {D n+1 - D n ) + 0(At 2 ) 


(3.4) 


where A, f?, and C are the Jacobian matrices 


■ dE , _dF . _aa 
'■* 3D' dD' dD 


(3.5) 


To evaluate these, recall that 
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Then 



where 


re* 

ZxP 

ZyP 

ZzP ' 

tx 

U + ( z u 

ZyU 

£z« 

Zy 

ZxV 

U+tyV 

ZzV 

Uz 

Zx™ 

Zy w 

U + fzttf J 


U = + ixU + ZyV + £ Z W 


(3.6) 


B and C are obtained in a similar manner, and all three can be represented by the following, 
where A t = A,B,ot C for i = 1,2, or 3, respectively. 


\L 0 

L i/3 

L 20 

UP ■ 

L 1 

Q + L\U 

L^u 

L 3 u 

L 2 

L\ v 

Q + L 2 v 

L 3 v 

L z 

L\w 

L 2 w 

Q + L 3 w . 


where 

Q — Lq + L\U ■+■ L>2 v "b L3W 

Lo = (£t)t> L>\ = (£i) xi L2 = (£i)y? B3 — (£t)z 

Substituting equation (3.4) into equation (3.3) results in the governing equation in delta form 


i + ^r+' 

2 


A tJ' 


^(i n - ro + 6„{B n - r 2 ) + s<{c n - r 3 )] | ( D n+l - D n ) 
6 t {E - E v ) n + 6r,(F - F v ) n + tf f (G - G v ) n 


(3.8) 


jn+ 1 

+ 1 -j- - > 1 


where 


r iD 


n-f- 1 


YoD 


n-fl _ £m+l 


T 3 D n+1 = G” +1 


At this point it should be noted that the notation of the form [^(A — T)]D refers to ^( AD ) — 
^ £ (r£>) and not § jD - § ^D. 


Approximate Factorization 

The solution of equation (3.8) would involve a formidable matrix inversion problem. 
With the use of an ADI (alternating direction implicit) type scheme, the problem could be reduced 
to the inversion of three matrices of small bandwidth, for which there exist some efficient solution 
algorithms. The particular ADI form used here is known as approximate factorization (ref. 24). 
Using this, the governing equation becomes 

L^L^D 71 * 1 - D n ) = RHS (3.9) 
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where 


L(= I+^J" +, 6 ( (A"-T,) 

L„= 1+ ^J n+l «,(B”-r,) (3.10) 

L, = I + ^J" +, 6 t (C" - r.) 

and RHS is the right-hand side of equation (3.8). When second-order central differencing is 
used, the solution to this problem becomes the inversion of three block tridiagonal matrices. The 
solution procedure is to solve these three problems 

(. Lr,)AD = RHS 

[L i )AD = AD (3.11) 

(LjA £> n+1 = A D. 

As an example of one of the inversions, the procedure is illustrated for the f-sweep, when (/ + 
\ ArJ n+1 6^C)AD = AD is solved for AD, and C = The computational domain indicies 

used are j, k, and 1 for the £, 77 , and f directions, respectively. Then using an LU decomposition, 
the following matrix equation is solved once for each j = 1 to jmax and k = 1 to kmax. 

I 0 

I 

0 -Cj tkt 


0 

0 

ADj,k,2 

A£>;, fc,3 

ADj ) k,lmax— 1 

0 

where 

C = -ArJ n+l (C n - T 3 ) 

4 

C is given by equation (3.7) with i = 3, and T 3 is given by equation (3.15). 

Using this implementation, there will be no change in the variables at the boundaries, 
and the boundary conditions can be implemented explicitly. It is possible, however, to implement 
the boundary conditions implicitly, and this will be discussed in a later section. 

The factorization has introduced the following second-order cross product term into the 

equation: 

h 2 \6^A6r,B + SnBS^C + 6 f CS e A] AD + Q(h 3 ) (3.13) 
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where 


a = a h - r u B = 5 n -r 2 , c = c n -r 3 , k=~j n+l 

To maintain the second order accuracy of the scheme, the added terms must be kept smaller 
than the original terms in the equations everywhere in the computational domain. This puts a 
restriction on the size of the artificial compressibility parameter /?. The proper choice of (3 will be 
discussed in a later section. In applying the approximate factorization scheme, it has been found 
that the stability of the scheme is dependent on the use of some higher-order smoothing terms. 
These are used to damp out higher frequency oscillations which arise in the solution because of 
the use of second-order central differencing, and will be discussed in the following section. 


Higher-Order Smoothing 

Higher-order smoothing terms are required to make the present algorithm stable. These 
terms help to damp out the higher order oscillations in the solution that are caused by the use of 
second-order central differencing. The smoothing term is derived in conjunction with an upwind 
finite-difference approximation to ^ which is given as 


for u < 0 


(3u -4u _ j+u _ 2 ) 

v „ = ; 2Ai for « > 0 

X ' (-3u ; +4u ?h1 -u ?+2 ) 

2Ax 

where uy = u(jAx). This equation can be written, for all u 

1 


(3.14) 


V x u 


+ ■ 


4Ax 
u 


[(3Uj - 4uy_i + U]-i) + ( — 3lly + 4uy + i - U J+2 )] 
[(3u ; - 4Uj_! + Uj _ 2 ) - ( — 3Uj + 4 u J + 1 - Uy + 2 )] 


(3.15) 


4|u| Ax 

Replacing the first [ j term of the above equation by a central difference formula results in 
1 


V T u 


2Ax 


(Uj + l _ u ]-l) + 


4|u| Ax 


(Uj -2 - 4u ; _i + 6uy - 4u J + i + Uj + 2 ). 


(3.16) 


This suggests that the proper form of the higher-order smoothing term should be that of the 
second term on the right-hand side of equation (3.16). This fourth-order term is represented 
using first-order upwind and downwind finite-difference operators 


(V i A i ) 2 u 


1 


4Ax 


(Uy_ 2 - 4uy _ j + 6 Uj - 4uy + j + Uy +2 ) 


(3.17) 


Similarly, by starting with a first-order upwind equation, a second-order smoothing term can be 
derived 

V x A x u = — (u } 1 — 2 Uj + tty — j ) (3.18) 

To compare how a small perturbation u' would act under the influence of each of these smoothing 
terms the following equation is used. Moving with the fluid, the perturbation under second-order 
smoothing is given by 

du' d 2 u' 
dt 6 dx 2 
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whose solution is 


u' = e~ caH (ae tax + be~ ,ax ) 
Under forth-order smoothing the perturbation is governed by 


du' dV 
dt 6 dx 4 


= 0 


whose solution is 

u' = e~ ta t (ae iax + be~ iax ) 

Under the second-order smoothing, the lower frequency perturbations (ct < 1) will be 
damped out faster than they will be with the fourth-order smoothing. For higher-frequency 
perturbations, the fourth-order smoothing will be more effective in damping the perturbation. 

To preserve the tridiagonal nature of the system, only second-order smoothing is used 
on the left-hand side of the equation, and fourth-order smoothing is used on the right-hand side. 
The equation to be solved now becomes 

I + hS^{A n — T i) + 

/ + «„(£» - r,) + <iV„A„] 
i + he, (c° - r 3 ) + (.VjA, ] (zr +1 - £>") 

= RHS - [(V £ A £ ) 2 + (V, A ,) 2 + (V t A t ) 2 ] D n 

where e, and c e are implicit and explicit smoothing parameters. There is a good discussion on 
the characteristics of these parameters by Chang et al (ref. 27). For this analysis they assumed 
a perturbation p' given by 

p'(f,r) = /(<V“ £ 

where f is the free-stream direction. It was found that in order to damp out the numerical 
fluctuations as time advances, the smoothing parameters must satisfy 


t, > e e 



(3.20) 


Employing discrete Fourier analysis, a necessary condition on these coefficients can be derived, 
namely 

e, > 2f, (3.21) 


However, the exact relation between these two coefficients can only be determined by a nonlinear 
stability analysis. For the results given in this paper, c, is taken to be 0.3 and e e is 0.1 for the 
smoothing on the velocity components. For the smoothing on the pressure, larg-r values are 
usually used. This will be discussed further in a following section. 


Boundary Conditions 

An important part of any computer code is the proper implementation of boundary 
conditions. The code must be capable of handling the several different types of boundaries 
encountered in numerical simulations which include solid surface, in-flow and out-flow, and far- 
field boundaries. 
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Solid Surface 


At a solid surface boundary, the usual no-slip condition is applied. In general the grid 
point adjacent to the surface will be sufficiently fine so that constant pressure normal to the 
surface in the viscous boundary layer can be assumed. For a f = constant surface, this can be 
expressed as 



(3.22) 


This boundary condition can be implemented either explicitly or implicitly. The explicit imple- 
mentation uses equation (3.12) for the matrix inversion and then the pressure on the boundary 
can be computed at the end of the time-step using a suitable finite-difference form of equation 
(3.22). The implicit implementation, however, will enhance the stability of the code. This can be 
done during the f sweep by replacing the first row of the matrix equation (3.12) with 


where 


IADj,k,i + cADj t k,2 = f 


•-1 

0 

0 

0- 


’PL = 2 - PL=1 ' 

0 

0 

0 

0 

A 

0 

0 

0 

0 

0 

/ = 

0 

. 0 

0 

0 

0. 


0 


(3.23) 


In-flow, Out-flow and Far-field Conditions 


The in-flow and out-flow boundary conditions for an internal flow problem and the far- 
field boundary conditions for an external flow problem can be handled in much the same way. 
The incoming flow for both problems can be set to some appropriate constant as dictated by the 
problem. For example, at the inlet to a pipe, the pressure can be set to a constant and the velocity 
profile can be specified to be uniform. The conditions at the downstream however, are the most 
difficult to provide. Simple upwind extrapolation is not well-posed. The best convergence rate is 
obtained if global mass is conserved. So to ensure the best results, the velocities and pressure are 
first updated using a second-order upwind extrapolation. For an exit at L = lmax this is 


Q 


n 

lmax 


1 y Azi J ^imai-2 

AZ 2 1 

A*! 1 


(3.24) 


where 


Az i — Zlmaz Zlmax — 1 
A^2 — %lma x “ %lmax — 2 


Then, these extrapolated velocities are integrated over the exit boundary to obtain the outlet 
mass flux 

rhout = f V n da. (3.25) 

J exit 

Then the velocity components are weighted by the mass flux ratio to conserve global mass 


y'n+l 



m out 


(3.26) 
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If nothing further is done to update the boundary pressure, this can lead to discontinuities in the 
pressure because momentum is not being conserved. A method of weighting the pressure by a 
momentum correction was presented by Chang et al (ref. 27). They use the following to obtain 
a pressure which corresponds to the mass weighted velocities 



(3.27) 


where W is the contravariant velocity given by Eq. (2.20). In obtaining this formula, it has 
been assumed that the streamlines near the exit plane are nearly straight. Any appreciable 
deviation will cause a discontinuity in the pressure and may lead to an instability. To avoid 
this, a momentum weighted pressure was used. This was obtained by integrating the momentum 
corrected pressure p n+1 and the extrapolated pressure p” across the exit 


/; +1 = f p n+l da 

J exit 

i; = f P h dd. 

J exit 

The final outlet pressure is then taken to be 



(3.28) 


(3.29) 


Using these downstream boundary conditions global conservation of mass and momentum are 
ensured and the scheme will not introduce instabilities into the flow field. 


PSEUDOCOMPRESSIBILITY 

In an incompressible flow, a disturbance in the pressure causes waves which travel with 
infinite speed. Introducing pseudocompressibility results in waves of finite speed, where the mag- 
nitude of the speed depends upon the pseudocompressibility constant ( 3 . In a true incompressible 
flow, the pressure field is affected instantaneously by a disturbance in the flow, but with pseu- 
docompressibility, there will be a time lag between the flow disturbance and its effect on the 
pressure field. In viscous flows, the behavior of the boundary layer is very sensitive to the stream- 
wise pressure gradient, especially when the boundary layer is separated. If separation is present, 
a pressure wave traveling with finite speed will cause a change in the local pressure gradient which 
will affect the location of the flow separation. This change in separated flow will feed back to the 
pressure field, possibly preventing convergence to a steady state. 


Characteristics of The Pressure Wave 

To analyze the behavior of this pressure wave and obtain the pseudo sound speed, a 
one-dimensional form of the governing equations is used 
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(4.1) 


dp du 

— + a — = o 

dt p dx 


du 
di + 


du 2 

dx 


= -^-r 

dx u 


where t w is the normalized shear stress at the wall. For the purpose of finding the pseudo sound 
speed, the shear stress is neglected. If the sound speed is denoted as c, and i denotes the position 
of the pressure wave, then an upstream traveling wave is moving at a speed of c — u. During a 
short time interval At, the wave will travel a short distance Ax. As At approaches zero, 



and 

du 

Yt = 

dp _ 

dt 

Substituting this into equation (4.1) and 


du dx 
dxdi = ' c 


«) 


du 

dx 

dp 


dp dx . . 

aim ={c ~ u) ai 


solving for c leads to 


c = YYTp 


The pseudo Mach number can be defined 


M 


v 

c 


u 

y/u 2 + /? 


(4.2) 


(4.3) 


(4.4) 


which is always less than one for any ft > 0. This is a necessary condition so that the pressure 
wave will propogate upstream and effect the whole flow field. 


Wave-Vorticity Interactions 


In their paper, Chang and Kwak (ref. 28) presented an analysis of the interaction 
between the finite-speed pressure wave and the development of the boundary layer in order to set 
a theoretical lower limit for j3. Their analysis is summarized here. The wave has to travel fast 
enough to allow proper distribution of the pressure while the boundary layer develops. In their 
analysis the following equations are obtained from equation (4.1) by linearizing the momentum 
flux locally, and cross differentiating 


which can be written as 


d 2 p - d 2 p d 2 p 

M* +2Uo di£i 

d 2 u - d 2 u d 2 tt 

~dt? +2U °dt8x ~ 


dx 

dr,,, 


dt 


at 


+ (U 0 + c) 


dx 



_ dr w I 
dt / 


(4.5) 


(4.6) 
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where c is the pseudo sound speed in equation (4.3) and Uq is the mean velocity in the primary 
flow direction. 

If the right-hand side of equation (4.6) is neglected, the result would be the characteristic 
equations for two families of linear waves 


d_ 

dt 


+ <u + c) az 


‘ d_ 

dt 


+ (u - c) 


d 

dx 



(4.7) 


The (+) family of waves propagate downstream with a velocity of c + u, and the ( ) family of 
waves propagate upstream with a velocity of c — u. 

If the shear stress term is included, then there is a coupling between the pseudo pressure 
waves and the vorticity spreading since the shear stress depends on the velocity field. Since the 
upstream traveling waves travel more slowly, these will be studied. For these waves, the coupling 
may be described by 


du 

~dt 


+ (u 


du dr 

= ~dt 


(4.8) 


To investigate the coupling effect, the problem of a flow through a channel of width x Tr j and 
length L is considered. The wave with the lowest wave number and therefore a length scale of L 
will be studied and so the following quantities will be used 


x = 


x 

r 


- = (u - cjt 



The rate of vorticity spreading for laminar flow is approximately given by 


d6 2 4 

W 

dt Re 

The scaled variable t„ was chosen such that 



(4.9) 


Using these variables, equation (4.8) becomes 


du du 
dt dx 


4L 


[i?e(u — c) J dt 


dr.. 


(4.10) 


where r„, = r,„(u, £(£„)). In this equation, the two u~ partial derivative terms are of order 1 and 
the shear stress partial derivative term is of the order of 1 or less, depending on tl.d size of the 
boundary layer. Therefore, the interactions between the waves and the vorticity can be decoupled 
if the coefficient on the right-hand side of equation (4.10) is very small 


4L 

Re(u — c) 


<< 1 


(4.11) 


By plugging the expression for the pseudo sound speed into this equation and assuming the 
normalized velocity of the primary flow u to be unity, the restriction for /? is obtained 
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(4.12) 


The physical meaning of this restriction can be thought of as follows. The time it takes an 
upstream wave to travel the distance of the computational domain of length L is 


Tr = 


t : - -it 


(4.13) 


Using equation (4.9), the time required for vorticity to spread a thickness 6 is approximatly 


Re ~ 
T„ « — S 2 

4 


(4.14) 


If it is required that the pressure wave travel the full length of the computational domain before 
the viscosity has diffused through the flow, then 


T l « T„ 


(4.15) 


By substituting in equations (4.13) and (4.14) into equation (4.15) and taking 6 — x rtJ = 1, the 
restriction in equation (4.11) is obtained. 

This analysis provides a lower bound for the choice of /?. For the upper bound of /?, the 
error introduced by the approximate factorization will be considered. Recalling equation (3.13), 
the form of the error is 


*=(f •/"+') [<«(*■- r,)«,(£“ - r) 


(4.16) 


This term must be kept smaller than the original terms in the equation. Including only the terms 
which contain /?, this restriction can be expressed as 


J' +l 6 e A" < 1 


(4.17) 


Recalling the expression for A" given by equation (3.7), the terms which have (3 in them give the 
following 


?».. m- 


(4.18) 


The term to the right of (3 in this inequality is essentially the change in -yr— in either the f,n, or c 

direction. An estimate of the order of magnitude of this term for the grids used in this study is 
given by 

°k(^)l“2 I*-' 9 ) 
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which puts the restriction on (3 


0(/3At) < 1 (4.20) 

For most problems, the restrictions for (3 given by equations (4.12) and (4.20) are satisfied with 
a value for / 3 in the range from 1 to 10. 

To test the validity of these restrictions, a test problem was run over a range of (3. The 
problem used was flow through a two-dimensional channel with a width of 1 and a length of 15. 
The prescribed inlet conditions were that of a uniform flow and a pressure of unity. The grid used, 
shown in figure 1, has 31 points in the cross-stream direction, and 65 points in the streamwise 
direction. The Reynolds number used is 1000 based on the duct width and inlet flow velocity; 
the time-step used in the calculations is At = 0.1. 

For a proper convergence, the above equations for the restrictions on (3 give the following 

range 

0.1236 < (3 < 10 (4.21) 

To test these values, the problem was run for the five values of (3 = 0.1, 1.0, 5.0, 10.0, and 50.0. 
The solution converged only for the cases of /? = 1.0, 5.0, and 10.0. The converged velocity profiles 
are shown for the case of (3 — 5.0 in figure 2. 

The convergence history of the five cases is shown in figure 3. The log (base 10) of the 
root-mean-square of the change in flow variables is plotted versus the computational time t . It 
can be seen that the calculations for the cases of /? = 0.1 and 50.0 become unstable within 100 
time steps and blow up, whereas the other cases converge to a stable solution. This indicates 
that the restrictions given in equations (4.12) and (4.20) are valid for this problem. In fact, these 
restrictions have proven to be valid for all cases yet encountered with this code. In figure 4, 
the effect of the various values of (3 on the incompressibility of the fluid is shown. The log of 
the root- mean-square of the divergence of the velocity field is plotted versus time. Again, the 
instability of the large and small (3 is seen, and the best convergence toward an incompressible 
solution is given with (3 = 5.0. 

To illustrate the effect of the traveling pressure wave in the channel flow, the following 
pressure contours are shown for the cases of (3 — 0.1 and /? - 5.0. In figure 5, the contours are 
shown for (3 = 5.0 at the computational times of 1.0, 2.5, and 4.0. The contours for )3 = 0.1 are 
shown at the same times in figure 6. In both cases, it can be seen that the pressure is dropping 
at the exit as the downstream boundary condition enforces conservation of global mass. For the 
j3 = 5.0 case, these contours propagate upstream as time progresses and the proper pressure 
drop begins to develop. For the (3 = 0.1 case, however, the pressure gradient does not propagate 
upstream nearly as fast. For this case at r = 2.5, the pressure gradient remains large near the 
outlet, the contours are not spreading upstream. The solution becomes unstable and this is shown 
by the interesting pattern in the contours as the pressure develops discontinuities. 


Effects of Smoothing on Pressure 

The explicit smoothing has a large effect on the convergence and accuracy of the pseudo- 
compressibility method. In particular, the explicit smoothing on the pressure can effect whether 
or not the solution converges to an incompressible flow field. In their paper, Chang and Kwak 
(ref. 28) showed that the pseudo pressure waves decay exponentially with time and vanish as the 
solution converges. Thus the change in pressure with time approaches zero. When there is no 
explicit smoothing added to the equation, as in equation (2.4), then the divergence of the velocity 
field identically approaches zero. However, when explicit smoothing is included, as the change in 
pressure approaches zero, the divergence of velocity approaches 

V-V - [(V £ A £ ) 2 + (V,A,) 2 + (V f A,) 2 ]p (4.22) 
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where 6,.(l) is the explicit smoothing parameter for the pressure. When the pressure gradients 
become substantial, as in the case when a region of separation is present, the smoothing term 
does not approach zero and this contaminates the incompressibility of the solution. 

To alleviate this problem, the value of 6, (1) is decreased with time. For the cases 
illustrated in the previous section, 6 C (1) was decreased by half of its value at every 50 time steps, 
from a starting value of 0.12, which is why there are kinks in the convergence curves in figures 3 
and 4. In this section, the same problem was used to illustrate the effect of the smoothing. For 
these cases, an exponentially decaying value of c, (1) was used 

e,(l) = 1.0e _, “ (4.23) 

for the values of a = 0.01, 0.05, 0.1, and 0.5. The results are presented in figures 7-10. Here, 
the logs of RMSDQ, RMSCO, and RMSDIV are plotted versus time, where RMSDQ is the 
root-mean-square of the change in the flow variables at every time step; where RMSCO is the 
root-mean-square of the right-hand side of the continuity equation, which includes the divergence 
of the velocity and the smoothing terms; and where RMSDIV is the root-mean-square of the 
divergence of the velocity field. 

Since the right-hand side of the continuity equation is equal to 

(-/JArdivV — smoothing terms) 

then when 0 = 5 and At = 0.1 and as the smoothing dies out, it is expected that RMSDIV 
will approach 2-RMSCO. In figure 7, for the case of a = 0.01, RMSCO rapidly converges but 
the divergence does not as it is being influenced heavily by the smoothing terms. In the next 
figures, it is seen that as the smoothing approaches zero, then the divergence of the velocity field 
converges at the same rate as RMSCO. 


COMPUTED RESULTS 

To test the accuracy of the code, two test problems are used. It was desired that both 
an internal flow and an external flow be tested. The first test case is that of laminar flow through 
a two-dimensional backward-facing step, the second is laminar flow around a two-dimensional 
circular cylinder. Both of these cases have experimental data available for comparison, and both 
have regions of separated flow which will test the ability of pseudocompressibility to distribute 
the pressure in the pressence of a shear layer without causing an instability. Since the cylinder 
flow is unstable for a Reynolds number greater than 40, the numerical solution will not converge 
and a psuedo unsteady behavior will be observed. 

Finally, a flow that is highly three-dimensional in nature is studied. For this, the flow 
about a three-dimensional cylinder on a flat plate is used. 


Backward Facing Step 

The flow over a backward-facing step was computed and the reattachment length was 
calculated for Reynolds numbers in the laminar range. The geometry of this problem is shown in 
figure 11. Here, the Reynolds number is based on twice the step height. The upstream boundary 
is located at the step and the inlet velocity profile there is prescribed to be a parabola (fully 
developed channel flow). The computations were conducted using a 65 x 33 grid similar to that 
shown in figure 1. The results of the experiment performed by Armaly et al (ref. 29) were used 
as a comparison to judge the accuracy of the computations. They reported that the flow at the 
step showed only a very small deviation from a parabolic profile, so that the placement of the 
inlet at the step should not introduce a significant error. 
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In the initial calculations, the code severely underpredicted the separation length at the 
higher range of Reynolds numbers. However, because of the presence of the explicit smoothing 
term in equation (4.22), the steady-state solution was not approaching a divergence free state. 
After the explicit smoothing parameter for pressure was decreased with time, the following results 
were obtained. In figure 12, the separation length versus the Reynolds number for three different 
computational codes and the experiment are plotted. 

INS3D shows remarkably good agreement with the experimental results for the full range 
of Reynolds numbers. The computations from the TEACH code were performed by Armaly et al 
(ref. 29) to coincide with their experimental results. The number of grid points they could use 
was limited by their computer’s core memory, and they claim this could explain the discrepancy 
between their numerical results and the experimental results for the higher Reynolds numbers. 
The data from Kim and Moin (ref. 21) are the results of running their fractional step code with 
the same boundary conditions and number of grid points that were used to run INS3D. Their 
code uses an approximate factorization scheme and solves a Poisson equation for pressure using 
a direct method based on trigonometric expansions. They use a staggered grid and the code is 
written for a uniform Cartesian mesh. 

In the INS3D calculations, opposite wall separation was present for Reynolds numbers 
of 500 and greater, which agrees with the results of the experiment. As in the case where the 
results were reported by the experiment, the opposite wall separation begins upstream from 
the reattachment point of the primary separation, and ends downstream from it. However, the 
magnitude of the computed opposite wall bubble is smaller than that found in the experiment. 
For the range of Reynolds numbers from 500 to 800, the length of this bubble was computed to be 
from 1 to 4 step-heights. The size in the experiment ranged from 5 to 8 step-heights in this range. 
Kim and Moin reported 7.8 step-heights for a Reynolds number of 600 and 11.5 step-heights for 
the 800 Reynolds number case. 

To compare the efficiency of the fractional step code with the INS3D code, the computing 
times used to solve the backward-facing step problem are presented in table 1. The convergence 
criteria were equivalent for the two cases. The computations were performed on a Cray X/MP 
and the times are given in units of CPU seconds. The calculation times for the two codes are 
comparable, with INS3D being a little more efficient. 

Table 1,'COMPUTING TIME FOR BACKSTEP PROBLEM. 


Reynolds Number 

INS3D 

Fractional Step 

100~ 

120 

92 

200 

121 

186 

300 

155 

315 

400 

230 

491 

500 

604 

649 

600 

614 

718 

700 

582 

803 

800 

601 

867 


In the next two figures, the results of INS3D’s calculations for a Reynolds r .mber of 700 
are shown. In figure 13, the streamlines are plotted, and in figure 14, velocity profiles at various 
points in the flow are shown. 
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Circular Cylinder 


The flow around a two-dimensional circular cylinder is rich in fluid physics. This flow 
is generally found to be steady for a Reynolds number less than 40 (ref. 30). Calculations for 
a Reynolds number of 40 have been carried out using INS3D and reported by Kwak et al (ref. 
25). The quantities they reported include the wake length, L„, 0 fc e ; the separation angle, 9 gep \ 
the pressure drag, Cd p ; and the pressure coefficients at the forward and rear stagnation points, 
Cp/and Cp r . All these quantities were in good agreement with a summary of experimental results 
found in literature and are repeated in table 2 for reference. 

Table 2r CIRCULAR CYLINDER FLOW AT RE = 40. 



Summary 

of 

literature 

INS3D 

C dp 

0.93-1.05 

1.03 

Cp/ 

1.14-1.23 

1.15 

C pr 

-0.47- -0.55 

-0.51 

v sep 

50°-53.9° 

52° 

Luiafce 

1.8-2. 5 

1.9 


Impulsively Started Symmetric Flow 

The development of the symmetric wake for an impulsively started circular cylinder was 
computed. The 120 by 85 grid shown in figure 15 was used for these calculations. The flow 
visualization results compare well with that of an experiment. In figure 16, the wake behind an 
impulsively started cylinder at a Reynolds number of 1200 is shown at nondimensional times of 
1.1, 1.3, 1.9, 2.4, 2.9, and 3.1 for both the experiment of Nagata et al (ref. 31) and the calculations 
of INS3D. The time-step used in the computations was 0.02. The similarity is striking, not only 
in the growth of the primary separation regions, but also in the appearances of two secondary 
bubbles near the point of separation. 

Bouard and Contanceau (ref. 32) also presented experimental results for an impulsively 
started circular cylinder. For the range 800 < Re < 5000, they noted what they call phenomenon 
a. First, for t > 1 , near the wall and about half-way between the stagnation and separation 
points a bulge appears forcing the streamlines out from the wall. They noted that this bulge 
develops into a little secondary eddy by t > 1.5. The very beginnings of this bulge can be seen 
in the computed streamlines at t = 1.3, and the secondary eddy is developed by t = 1.9. In 
the experiment they noted that the eddy grows to touch the outer boundary of the primary 
eddy, splitting off another small eddy from the primary eddy. This new eddy is equal in size and 
strength to the secondary eddy, thus there is a pair of secondary eddies formed. This can clearly 
be seen in the next three pictures of the computed streamlines. These same features have also 
been experimentally observed by Honji & Taneda (ref. 33), and Taneda (ref. 34), and numerically 
by Thoman & Szewczyk (ref. 35), and more recently by Loc (ref. 36), and Lecointe & Piquet 
(ref. 37). 

Vortex Shedding 

Even though the present algorithm is designed for obtaining steady-state solutions, 
pseudo-unsteady characteristics are examined in this section. In order to observe periodic behav- 
ior in the wake of the cylinder, it was necessary to introduce a perturbation in the flow field. For 
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this, surface roughness equal to 2 % of the diameter was applied to one side of the body. Once 
the periodic motion had started, no further disturbance was needed to continue to shedding, and 
a smooth grid was used for further calculations. The calculations were carried out at a Reynolds 
number of 1200, for which a Strouhal number in the range of 0.17 to 0.2 (ref. 30) has been found 
in experiment. Using a time-step of 0.01, the calculated Strouhal number varied from 0.357 to 
0.303 depending on what mesh was used. The value of 0.303 was obtained on the finest mesh 
used, that being 140 by 120. The lift and coefficient of pressure drag versus time are shown in 
figures 17 and 18, respectively, starting at a time when the smooth grid was used. These are the 
calculations for the 120 by 85 grid and the Strouhal number based on the frequency of the pres- 
sure drag over this time interval is 0.357. Since the calculations never approach a steady state, 
the pseudo pressure waves never die out, and the calculated flow will not become completely 
divergence-free. For these calculations, the root mean square of the divergence of the velocity 
field had an average value of 0.02, so it is not surprising that the calculations did not predict the 
correct Strouhal number. However, the phenomenon of vortex shedding is qualitatively shown by 
this method of computation. 

Streamlines at different times in the period are shown in figure 19 at times of 10.4, 11.1, 
11.8, 12.5, and 13.3. At r = 10.4, both the lift and drag are at a maximum, and it can be seen 
that a vortex is about to be shed from the upper half of the cylinder. Here it can be seen that 
the core streamlines on the top surface fully encircle the separated region and travel all the way 
to the lower half before starting downstream. At the next time frame, the separated region has 
begun to open up, and the lift starts to drop but it is still positive. In the next frame, the lift is 
nearly zero, and the drag is at a minimum. In the next frame the lower separated region begins 
to form and the lift is becoming negative, while the drag begins to increase. In the last frame, 
the lift has just turned past the minimum, and the streamlines are almost a reversed image of 
the first frame. 


Three-Dimensional Calculations 

The code INS3D has been used to calculate the flow for several different complex three- 
dimensional objects including a geometry used to model the Space Shuttle Main Engine (refs 25 
and 27). As an example of the three-dimensional capability, the flow is calculated for the problem 
of two parallel plates with a cylinder normal to the plate surface attached between the two 
plates. This geometry is similar to a flow in the Space Shuttle Main Engine where hot hydrogen 
gas under high pressure is flowing around posts which carry liquid oxygen. Since the flow is 
symmetric about the center plane between the plates, only one plate with half the cylinder will 
be used in the calculations. The inflow is prescribed to be a parabolic profile, and the Reynolds 
number is based on the average inflow velocity and the diameter of the cylinder. The grid and 
inflow are shown in figure 20, with a close up of the grid spacing at the cylinder-plate junction 
shown in figure 21. 

The solution to this problem was found to be steady, and some interesting features 
of the flow were found in the calculations. A saddle point of separation and a nodal point of 
reattachment appeared in front of the cylinder as shown in figures 22 and 23. Figure 22 is a top 
view of the velocity vectors in the z plane next to the wall at z/D = 0.02, where D is the diameter 
of the cylinder. The saddle point of separation appears at y/D = -2.9, and the nodal point of 
reattachment is at y/D = -0.56. Figure 23 is a side view of the velocity vectors in the plane at 
x/D - 0.0. Figure 24 shows particle traces of the computed flow for particles released on either 
side of the x/D = 0.0 plane. Particles which are released near the region of the separated flow 
between the saddle point of separation and the nodal point of reattachment show the presence of 
a horseshoe vortex. Particles released between the nodal point of separation and the body show 
the existance of a pair of spiraling vortex filaments which wrap around the cylinder and spiral 
upward. This secondary vortex has a sense opposite that of the horseshoe vortex. 
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CONCLUSION 

An incompressible flow solver has been developed which is capable of computing the 
steady-state solution to a realistic three-dimensional geometry. Pseudocompressibility has been 
found to be an efficient way to solve the pressure field, and although this limits the time-accuracy, 
it is capable of showing the qualitative nature of an unsteady problem. The ability of this method 
to perform well depends heavily on the proper choice of the pseudocompressibility constant 0 . 
This paper has presented an analysis of the effect on the calculations of using different values 
of 0 , and on this basis has presented guidelines for the proper choice of 0 . These guidelines are 
substantiated with example calculations. The accuracy of the method is affected by the amount of 
smoothing used, and whereas the algorithm requires some smoothing for stability, by decreasing 
the smoothing as a steady state is reached, accurate solutions to the incompressible Navier-Stokes 
equations can be obtained. The computing time has been shown to be resonable when compared 
to another method which used a similar algorthim for its solution process. Finally, the code was 
shown to be able to solve a three-dimensional problem with some very interesting and complex 
fluid phenomena. 

Such a code could be a valuable tool for many enginering problems, particularly if some 
more study of the time-accuracy of the method were done. It is possible that an iterative-type 
scheme could help reduce the divergence of the velocity field at the end of each time-step. Such 
a scheme should be similar to the fractional step method iterations shown in equations (1.1) and 
(1.2). This iteration would be very much like artificial compressibility in that there would be 
artificial pressure waves distributing the pressure. If the iterations were carried out until they 
converged, this would correspond to the waves dying out completely, and the flow would be 
divergence-free. However, instead of iterating until the flow field is completely divergent-free, the 
iterations could be allowed to stop after the divergence of the velocity has been reduced by an 
order of magnitude or so. Through testing, a criteria could be found as to how much iterating 
must be done to get a desired level of quantitative accuracy. In other words, the way for the 
artifical compressibility to obtain time-accurate answers is to let the pressure waves travel further 
for each time-step than they do in the current one-step formulation. 
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Figure 1 Grid used for channel flow. 



Figure 2 Velocity profiles for /? = 5.0. 
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Figure 7 Effect of decreasing smoothing with time a — 0.01. 



Figure 8 Effect of decreasing smoothing with time a = 0.05. 
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Figure 9 Effect of decreasing smoothing with time a = 0.1. 



Figure 10 Effect of decreasing smoothing with time a = 0.5. 
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Figure 11 Flow over a backward-facing step. 
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Figure 12 Reattachment length versus Reynolds number. 
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Figure 13 Streamlines for backstep flow at Re = 700. 



0 5 10 15 20 

X 

Figure 14 Velocity vectors for backstep flow at Re = 700. 
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Figure 15 Grid for a circular cylinder. 
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T=1.9 


Figure 16a Experimental pictures and computational streamlines for an impulsively started cylinder 
at r -- 1.1, 1.3, and 1.9. 
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Figure 16b Experimental pictures and computational streamlines for an impulsively started cylinder 
at t = 2.4, 2.9, and 3.1. 
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Figure 17 Lift versus time for periodic flow over a cylinder. 
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Figure 18 Coefficient of pressure drag versus time. 


38 


I 



Figure 19 Streamlines of periodic flow past a cylinder. 
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Figure 21 Grid spacing near cylinder-plate junction. 
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